library(gsynth)
library(panelView)
library(DataCombine)
dat <- read.csv("prepped-gensynth.csv")

#### Replicates Table 2 and Figure 1 in the main paper. Also does some robustness
#### checks and other visualizations. Note that exact numbers will vary from run
#### to run becasue of random variation in the bootstrap.

#### NOTEs:
## We used gensynth 1.0.9
##
## install.packages("devtools")
## library(devtools)
## install_github("xuyiqing/gsynth")
##
## You need to set the SDs variable (lines 22, 93) in order to look at different shock 
## definitions. Set for paper figures at the moment.

#### Define intervention

## Take a lookuniq
summary(dat$v2x_polyarchy_diff[dat$v2x_polyarchy_diff > 0])
summary(dat$v2x_polyarchy_diff[dat$v2x_polyarchy_diff > 0.05])
summary(dat$v2x_polyarchy_diff[dat$v2x_polyarchy_diff > 0.1])

## mean + 2sd jump in diff
SDs <- 2
THRESH <- mean(dat$v2x_polyarchy_diff, na.rm=T) + SDs*sd(dat$v2x_polyarchy_diff, na.rm=T)
ccs <-  unique(dat$countrycode)
hits <- sapply(ccs, function (cc)
  sapply(dat$v2x_polyarchy_diff[dat$countrycode==cc], function (ts)
    ts > THRESH))
hits2<- lapply(hits, function (x) {
  cur <- FALSE
  sapply(x, function (y) {
    if (is.na(y)) y <- FALSE
    if (y) cur <<- TRUE
    cur
  })
})

dat$treated <- dat$v2x_polyarchy_diff > THRESH
for (i in 1:length(hits2)) dat$treated[dat$countrycode==ccs[i]] <- hits2[[i]]

## Visualize
dat$treated <- as.numeric(dat$treated)
panelView(AID100 ~ treated, data=dat, index=c("countrycode", "year"))
panelView(AID100 ~ treated, data=dat, index=c("countrycode", "year"), 
          type='outcome', ylim=c(0, 100))

### Fit model with gov aid
dat <- slide(dat, Var="AID100", GroupVar="countrycode", slideBy=-1)
dat$AID100.l1 <- dat$"AID100-1"

fit <- gsynth(AID100 ~ treated + PRF01 + DIF07 + FPP01 + POL05 + POL25 + AID100.l1,
              min.T0=9,
              data=dat[!is.na(dat$PRF01) & !is.na(dat$AID100.l1),],
              index=c("countrycode", "year"),
              force="two-way", CV=TRUE, r=c(0,7),
              inference="parametric", nboots=100, se=TRUE)

## This is Figure 1, panel B
plot(fit) 

plot(fit, type = "counterfactual", raw = "all", main="")
plot(fit, type="raw", ylim=c(0,10))

## No LDV for reference
fit <- gsynth(AID100 ~ treated + PRF01 + DIF07 + FPP01 + POL05 + POL25,
              min.T0=9,
              data=dat[!is.na(dat$PRF01),],
              index=c("countrycode", "year"),
              force="two-way", CV=TRUE, r=c(0,7),
              inference="parametric", nboots=100, se=TRUE)

### Now from aid to polyarchy.

## Reset data
rm(list=ls())
dat <- read.csv("../data/prepped-gensynth.csv")

dat <- slide(dat, Var="AID100", GroupVar="countrycode", slideBy=-1)
dat$AID100.l1 <- dat$"AID100-1"
dat$AID100_diff <- dat$AID100 - dat$AID100.l1

SDs <- 1
THRESH <- mean(dat$AID100_diff, na.rm=T) + SDs * sd(dat$AID100_diff, na.rm=T)
ccs <-  unique(dat$countrycode)
hits <- sapply(ccs, function (cc)
  sapply(dat$AID100_diff[dat$countrycode==cc], function (ts)
    ts > THRESH))
hits2<- lapply(hits, function (x) {
  cur <- FALSE
  sapply(x, function (y) {
    if (is.na(y)) y <- FALSE
    if (y) cur <<- TRUE
    cur
  })
})

dat$treated <- dat$AID100_diff > THRESH
for (i in 1:length(hits2)) dat$treated[dat$countrycode==ccs[i]] <- hits2[[i]]

dat <- slide(dat, Var="v2x_polyarchy", GroupVar="countrycode", slideBy=-1)
dat$v2x_polyarchy.l1 <- dat$"v2x_polyarchy-1"

panelView(v2x_polyarchy ~ treated, data=dat, index=c("countrycode", "year"))

fit2 <- gsynth(v2x_polyarchy ~ treated + PRF01 + DIF07 + FPP01 + POL05 + POL25 
               + v2x_polyarchy.l1,
               min.T0=9,
               data=dat[!is.na(dat$PRF01) & !is.na(dat$v2x_polyarchy.l1),],
               index=c("countrycode", "year"),
               force="two-way", CV=TRUE, r=c(0,7),
               inference="parametric", nboots=100, se=TRUE)

## This is Figure 1, panel A
plot(fit2)

plot(fit2, type = "counterfactual", raw = "none", main="")

## No LDV version
fit2 <- gsynth(v2x_polyarchy ~ treated + PRF01 + DIF07 + FPP01 + POL05 + POL25,
               min.T0=9,
               data=dat[!is.na(dat$PRF01) & !is.na(dat$v2x_polyarchy),],
               index=c("countrycode", "year"),
               force="two-way", CV=TRUE, r=c(0,7),
               inference="parametric", nboots=100, se=TRUE)
